set scheme sj
// double threshold model
sjlog using xthreg1,replace
use hansen1999, clear
*tabstat i q1 c1 d1, stat(min p25 p50 p75 max) col(stat)
* single threshold model
xthreg i q1 q2 q3 d1 qd1, rx(c1) qx(d1) thnum(1) grid(400) trim(0.01) bs(300) 
sjlog close,replace

* triple threshold model given the result above
sjlog using xthreg2, replace
xthreg i q1 q2 q3 d1 qd1, rx(c1) qx(d1) thnum(3) grid(400) ///
trim(0.01 0.01 0.05) bs(0 300 300) thgiven nobslog noreg
sjlog close, replace

* plot the LR statistic vs threshold
sjlog using xthreg3, replace
_matplot e(LR21), columns(1 2) yline(7.35, lpattern(dash)) ///
	connect(direct) msize(small) mlabp(0) mlabs(zero) ///
	ytitle("LR Statistics") xtitle("First Threshold") recast(line) name(LR21) nodraw
_matplot e(LR22), columns(1 2) yline(7.35, lpattern(dash)) ///
	connect(direct) msize(small) mlabp(0) mlabs(zero) ///
	ytitle("LR Statistics") xtitle("Second Threshold") recast(line) name(LR22) nodraw
graph combine LR21 LR22, cols(1)	
sjlog close, replace

// estimate triple thrshold model directly
sjlog using xthreg4, replace
xthreg i q1 q2 q3 qd1 d1, rx(c1) qx(d1) thnum(3) grid(400) ///
trim(0.01 0.01 0.05) bs(300 300 300) 
sjlog close,replace
